

* First set your Stata local working directory to where you unpacked the repository 
* use EASmqsREVISION.dta, clear

* save on some memory
drop if hchamber==1
drop if hnum==71 | hnum==78 | hnum==84 | hnum==85
drop if hmedicaid==1

gen speccte=0
replace speccte=1 if cnum==101 | cnum==106 | cnum==109 | cnum==111
* 101==W&M/H_107th, 106==E&C/H_107th, 109==C/H&E_106th, 111==W&M/H_106th -- these are special committees

*create count dependent variables assuming selection - if the member doesn't as any questions at the hearing then assume not present
gen nauth_sel=nauth
replace nauth_sel=. if totalmqs==0
gen nanecdist_sel=nanecdist
replace nanecdist_sel=. if totalmqs==0

* reshape data - needed for stacking the DVs
gen vartype1=1  /* participate */
gen vartype2=2  /* experiential question count */
gen vartype3=3  /* analytical question count */
gen vartype4=4  /* lnhltmoney */
gen vartype5=5  /* lnallmoney */
gen vartype6=6  /* xpctemphmean */
gen vartype7=7  /* resorgct */
ren obs id
reshape long vartype, i(id) 

* create dummy variables to indicate equation rows
gen o1=cond(vartype==1,1,0)
gen o2=cond(vartype==2,1,0)
gen o3=cond(vartype==3,1,0)
gen h1=cond(vartype==4,1,0)
gen h2=cond(vartype==5,1,0)
gen w1=cond(vartype==6,1,0)
gen w2=cond(vartype==7,1,0)

* NOTE: "x" prefix means the variable was standardized; "ln" in the variable name indicates the variable was logged

* create stacked dependent variable
gen resp=participate if o1==1
replace resp=nanecdist_sel if o2==1
replace resp=nauth_sel if o3==1
replace resp=xlnhltmoney if h1==1
replace resp=xlnallmoney if h2==1
replace resp=xpctemphmean if w1==1
replace resp=resorgct if w2==1


* put independent variables in appropriate rows. See paper for variable definitions. Not all of these are used in the model.
gen xlnhltmoney_o1=xlnhltmoney*o1
gen xlnallmoney_o1=xlnallmoney*o1
gen chair_o1=chair*o1
gen republican_o1=republican*o1
gen xsmedicare_o1=xsmedicare*o1
gen xyearsofservice_o1=xyearsofservice*o1
gen xmostrecent_o1=xmostrecentmargin*o1
gen xdwnomf_o1=xdwnomf*o1
gen hltoccup_o1=hltoccup*o1
gen xoasdibenx_o1 =xoasdibenx*o1
gen xindemphltx_o1 =xindemphltx*o1
gen xcqmdnincm_o1=xcqmdnincm*o1
gen xnewsitems_o1=xnewsitems*o1
gen hmedicaid_o1=hmedicaid*o1
gen year2000_o1=year2000*o1
gen xlength_o1=xlength*o1
gen xf1expertarg_o1=xf1expertarg*o1    /* for possible use for starting values.... */
gen xresorgct_o1=xresorgct*o1
gen xpctemphmean_o1=xpctemphmean*o1
gen speccte_o1=speccte*o1

gen xlnhltmoney_o2=xlnhltmoney*o2
gen xlnallmoney_o2=xlnallmoney*o2
gen chair_o2=chair*o2
gen republican_o2=republican*o2
gen xsmedicare_o2=xsmedicare*o2
gen xyearsofservice_o2=xyearsofservice*o2
gen xmostrecent_o2=xmostrecentmargin*o2
gen xdwnomf_o2=xdwnomf*o2
gen hltoccup_o2=hltoccup*o2
gen xoasdibenx_o2 =xoasdibenx*o2
gen xindemphltx_o2 =xindemphltx*o2
gen xcqmdnincm_o2=xcqmdnincm*o2
gen xnewsitems_o2=xnewsitems*o2
gen hmedicaid_o2=hmedicaid*o2
gen year2000_o2=year2000*o2
gen xlength_o2=xlength*o2
gen xf1expertarg_o2=xf1expertarg*o2
gen xresorgct_o2=xresorgct*o2
gen xpctemphmean_o2=xpctemphmean*o2
gen speccte_o2=speccte*o2

gen xlnhltmoney_o3=xlnhltmoney*o3
gen xlnallmoney_o3=xlnallmoney*o3
gen chair_o3=chair*o3
gen republican_o3=republican*o3
gen xsmedicare_o3=xsmedicare*o3
gen xyearsofservice_o3=xyearsofservice*o3
gen xmostrecent_o3=xmostrecentmargin*o3
gen xdwnomf_o3=xdwnomf*o3
gen hltoccup_o3=hltoccup*o3
gen xoasdibenx_o3=xoasdibenx*o3
gen xindemphltx_o3=xindemphltx*o3
gen xcqmdnincm_o3=xcqmdnincm*o3
gen xnewsitems_o3=xnewsitems*o3
gen hmedicaid_o3=hmedicaid*o3
gen year2000_o3=year2000*o3
gen xlength_o3=xlength*o3
gen xf1expertarg_o3=xf1expertarg*o3
gen xresorgct_o3=xresorgct*o3
gen xpctemphmean_o3=xpctemphmean*o3
gen speccte_o3=speccte*o3

gen xctehltmnoi_h1=xctehltmnoi*h1
gen chair_h1=chair*h1
gen republican_h1=republican*h1
gen xsmedicare_h1=xsmedicare*h1
gen xyearsofservice_h1=xyearsofservice*h1
gen xmostrecent_h1=xmostrecentmargin*h1
gen xdwnomf_h1=xdwnomf*h1
gen hltoccup_h1=hltoccup*h1
gen xoasdibenx_h1=xoasdibenx*h1
gen xindemphltx_h1=xindemphltx*h1
gen xcqmdnincm_h1=xcqmdnincm*h1
gen speccte_h1=speccte*h1

gen xcteallmnoi_h2=xcteallmnoi*h2
gen chair_h2=chair*h2
gen republican_h2=republican*h2
gen xsall_h2=xsall*h2
gen xyearsofservice_h2=xyearsofservice*h2
gen xmostrecent_h2=xmostrecentmargin*h2
gen xdwnomf_h2=xdwnomf*h2
gen hltoccup_h2=hltoccup*h2
gen xoasdibenx_h2=xoasdibenx*h2
gen xindemphltx_h2=xindemphltx*h2
gen xcqmdnincm_h2=xcqmdnincm*h2
gen speccte_h2=speccte*h2

gen xhicost_w1=xhicost*w1
gen xlnhltmoney_w1=xlnhltmoney*w1
gen chair_w1=chair*w1
gen republican_w1=republican*w1
gen xsmedicare_w1=xsmedicare*w1
gen xyearsofservice_w1=xyearsofservice*w1
gen xmostrecent_w1=xmostrecentmargin*w1
gen xdwnomf_w1=xdwnomf*w1
gen hltoccup_w1=hltoccup*w1
gen xoasdibenx_w1=xoasdibenx*w1
gen xindemphltx_w1=xindemphltx*w1
gen xcqmdnincm_w1=xcqmdnincm*w1
gen xnewsitems_w1=xnewsitems*w1
gen hmedicaid_w1=hmedicaid*w1
gen year2000_w1=year2000*w1
gen xlength_w1=xlength*w1
gen xctehltmnoi_w1=xctehltmnoi*w1
gen xresorgct_w1=xresorgct*w1
gen speccte_w1=speccte*w1

gen xhicost_w2=xhicost*w2
gen xlnhltmoney_w2=xlnhltmoney*w2
gen chair_w2=chair*w2
gen republican_w2=republican*w2
gen xsmedicare_w2=xsmedicare*w2
gen xyearsofservice_w2=xyearsofservice*w2
gen xmostrecent_w2=xmostrecentmargin*w2
gen xdwnomf_w2=xdwnomf*w2
gen hltoccup_w2=hltoccup*w2
gen xoasdibenx_w2=xoasdibenx*w2
gen xindemphltx_w2=xindemphltx*w2
gen xcqmdnincm_w2=xcqmdnincm*w2
gen xnewsitems_w2=xnewsitems*w2
gen hmedicaid_w2=hmedicaid*w2
gen year2000_w2=year2000*w2
gen xlength_w2=xlength*w2
gen xctehltmnoi_w2=xctehltmnoi*w2
gen speccte_w2=speccte*w2

* these assign latent variables to equations, the first for eta1 and the second for eta2 (e.g., eta2 not in w1 and w2 equations etc.)
eq fo3o2o1w1w2: o3 o2 o1 w1 w2
eq fo3o2o1h1h2: o3 o2 o1 h1 h2

* these starting values speed things up by feeding in the solution.
mat startv=[-.06826186, .38728949, .09680224, -.18501639, .04645411, .07696352, -.14961891, .08054587, .10770327, .0835075, .16391579, .1322913, -.36697578, -.54619201, -.13636052, .0610982, .17639153, -.02690227, -.07089303, .29085928, .28447353, .2362559, -.17605898, .05368066, .61736794, -2.5213567, -1.0574422, -.03594962, .84270534, .17761401, .50240525, .55691917, .20571062, .21324267, -.19132674, -.13119932, -.25047695, -.05570252, .14419407, .30267289, -.2116397, -.37895974, .11725674, .28935496, .13282184, -.17771969, .0850025, -.03481522, 1.5466113, -.80470301, -.1507387, -.69154355, -1.2659989, -1.1124714, .0273119, .88663009, 2.0537585, 1.0830145, .24733337, .77791715, .05981288, -.60305886, 2.4310575, -.83585808, -1.0737716, 1.934085, -.40708813, .67431231, .06670512, 1.0869563, -2.455722, 1.0983795, 1.6570756, 1.3790062, .44428856]

* for how to use GLLAMM, see https://biostat.jhsph.edu/~fdominic/teaching/bio656/software/gllamm.manual.pdf pp. 20-27
* you can ignore the overflow/concavity warning messages -- concavity was not an issue for estimation on the Stata version used for the manuscript and the results in Stata 18 replicate the manuscript results exactly.
gllamm resp xresorgct_o1 xresorgct_o2 xresorgct_o3 xmostrecent_o1 xmostrecent_o2 xmostrecent_o3 xyearsofservice_o1 xyearsofservice_o2 xyearsofservice_o3 xdwnomf_o1 xdwnomf_o2 xdwnomf_o3 republican_o1 republican_o2 republican_o3 xoasdibenx_o1 xoasdibenx_o2 xoasdibenx_o3 xnewsitems_o1 xnewsitems_o2 xnewsitems_o3 xsmedicare_o1 xsmedicare_o2 xsmedicare_o3 xpctemphmean_o1 xpctemphmean_o2 xpctemphmean_o3 xlength_o1 xlength_o2 xlength_o3 xctehltmnoi_h1 chair_h1 republican_h1 xsmedicare_h1 xmostrecent_h1 h1 xcteallmnoi_h2 chair_h2 republican_h2 xsall_h2 xmostrecent_h2 h2 speccte_w1 xctehltmnoi_w1 xnewsitems_w1 year2000_w1 xlength_w1 w1 speccte_w2 xctehltmnoi_w2 xnewsitems_w2 year2000_w2 xlength_w2 w2 xlnhltmoney_o3 chair_o3 o3 xlnhltmoney_o2 chair_o2 o2 xlnhltmoney_o1 xlnallmoney_o1 chair_o1 o1, i(id mnum) eq(fo3o2o1w1w2 fo3o2o1h1h2) nrf(1 1) family(binom poisson poisson gauss gauss gauss poisson) link(probit log log ident ident ident log) fv(vartype) lv(vartype) nocons nip(8) trace adapt from(startv) copy


* optionally, save estimated results
* estsave, gen(model_8ip) label(results for APSR 2007)
* keep model_8ip
* save EASmqsREVISION_results.dta, replace

* you can load the saved results by uncommenting these lines:
* use EASmqsREVISION_results.dta, clear
* install estsave if you haven't already, and then to see the results open the results file and run:
* estsave, from( model_8ip )



*generate three graphs of results:

* final graph latent variable posterior estimates -- you can read about gllapred in  https://biostat.jhsph.edu/~fdominic/teaching/bio656/software/gllamm.manual.pdf pp. 27-29
capture drop _merge
gllapred fac1, u
bysort mnum: gen mkeeper1=(_n==1)
drop if mnum==.
gsort +fac1m2 -mkeeper1
gen rank1=sum(mkeeper1)
gen upper1=fac1m2+1.96*fac1s2+.1
serrbar fac1m2 fac1s2 rank1 if mkeeper1==1, scale(1.96) xtitle("Analytical Capacity Rank") ytitle("Analytical Capacity Expected Value") mvopts(mcolor(black) msize(small)) legend(off) yline(0, lwidth(thin)) lcolor(black) lwidth(thin) scheme(s1manual) 

* final graph analytical question count and analytical capacity
gen aqc=exp(2.054-.136+fac1m2)
gen aqc_upper=exp(2.054-.136+fac1m2+1.96*fac1s2)
gen aqc_lower=exp(2.054-.136+fac1m2-1.96*fac1s2)
gen eqc=exp(.778-.546+(-2.456*fac1m2))
gen eqc_upper=exp(.778-.546+(-2.456*(fac1m2+1.96*fac1s2)))
gen eqc_lower=exp(.778-.546+(-2.456*(fac1m2-1.96*fac1s2)))
twoway lowess aqc_lower fac1m2, sort lpattern(dot) lwidth(medium) || lowess aqc_upper fac1m2, sort lpattern(dot) lwidth(medium) || lowess aqc fac1m2, sort lwidth(medium) lcolor(black) || lowess eqc_lower fac1m2, sort lpattern(shortdash) lwidth(thin) || lowess eqc_upper fac1m2, sort lpattern(shortdash) lwidth(thin) || lowess eqc fac1m2, sort lwidth(medium) lcolor(black) ||, plotregion(style(none)) legend(off) xtitle("Analytical Capacity") ytitle("Expected Question Count") scheme(s1manual)

*final graph contributions 
gen hc=-.131+.206+(1.657*fac1m2)
gen hc_upper=-.131+.206+(1.657*(fac1m2+1.96*fac1s2))
gen hc_lower=-.131+.206+(1.657*(fac1m2-1.96*fac1s2))
gen ac=-.379+.144+(1.379*fac1m2)
gen ac_upper=-.379+.144+(1.379*(fac1m2+1.96*fac1s2))
gen ac_lower=-.379+.144+(1.379*(fac1m2-1.96*fac1s2))
gen ehc=exp(hc*1.033+11.13)/1000     /* converts xlnhltmoney back into 1000s dollars -- lnhltmoney mean=11.13, sd=1.033 */
gen ehc_upper=exp(hc_upper*1.033+11.13)/1000
gen ehc_lower=exp(hc_lower*1.033+11.13)/1000
gen eac=exp(ac*.747+13.46)/1000     /* converts xlnallmoney back into 1000s dollars -- lnallmoney mean=13.46 sd=0.747 */
gen eac_upper=exp(ac_upper*.747+13.46)/1000
gen eac_lower=exp(ac_lower*.747+13.46)/1000
twoway lowess ehc_lower fac1m2, sort lpattern(dot) lwidth(medium) || lowess ehc_upper fac1m2, sort lpattern(dot) lwidth(medium) || lowess ehc fac1m2, sort lwidth(medium) lcolor(black)  || lowess eac_lower fac1m2, sort lpattern(shortdash) lwidth(thin) || lowess eac_upper fac1m2, sort lpattern(shortdash) lwidth(thin) || lowess eac fac1m2, sort lwidth(medium) lcolor(black) ||, plotregion(style(none)) legend(off) xtitle("Analytical Capacity") ytitle("Expected Contributions (in $1000s)") scheme(s1manual)



/*
e(b)[1,75]
                  resp:               resp:               resp:               resp:               resp:               resp:
          xresorgct_o1        xresorgct_o2        xresorgct_o3      xmostrecent_o1      xmostrecent_o2      xmostrecent_o3
y1          -.06826209           .38728948           .09680222          -.18501686           .04645424           .07696354

                  resp:               resp:               resp:               resp:               resp:               resp:
    xyearsofservice_o1  xyearsofservice_o2  xyearsofservice_o3          xdwnomf_o1          xdwnomf_o2          xdwnomf_o3
y1           -.1496194           .08054638            .1077034           .08350773           .16391531           .13229114

                  resp:               resp:               resp:               resp:               resp:               resp:
         republican_o1       republican_o2       republican_o3       xoasdibenx_o1       xoasdibenx_o2       xoasdibenx_o3
y1          -.36697611          -.54619127          -.13636021           .06109834           .17639156           -.0269023

                  resp:               resp:               resp:               resp:               resp:               resp:
         xnewsitems_o1       xnewsitems_o2       xnewsitems_o3       xsmedicare_o1       xsmedicare_o2       xsmedicare_o3
y1          -.07089367           .29085933           .28447351           .23625676          -.17605952           .05368054

                  resp:               resp:               resp:               resp:               resp:               resp:
       xpctemphmean_o1     xpctemphmean_o2     xpctemphmean_o3          xlength_o1          xlength_o2          xlength_o3
y1            .6173727          -2.5213571           -1.057442          -.03595011           .84270529            .1776139

                  resp:               resp:               resp:               resp:               resp:               resp:
        xctehltmnoi_h1            chair_h1       republican_h1       xsmedicare_h1      xmostrecent_h1                  h1
y1           .50240522           .55691928           .20571057           .21324269          -.19132675          -.13119932

                  resp:               resp:               resp:               resp:               resp:               resp:
        xcteallmnoi_h2            chair_h2       republican_h2            xsall_h2      xmostrecent_h2                  h2
y1            -.250477          -.05570241           .14419402           .30267292          -.21163971          -.37895974

                  resp:               resp:               resp:               resp:               resp:               resp:
            speccte_w1      xctehltmnoi_w1       xnewsitems_w1         year2000_w1          xlength_w1                  w1
y1           .11725655           .28935513           .13282188          -.17771939           .08500248          -.03481523

                  resp:               resp:               resp:               resp:               resp:               resp:
            speccte_w2      xctehltmnoi_w2       xnewsitems_w2         year2000_w2          xlength_w2                  w2
y1           1.5466102           -.8047022          -.15073867          -.69154351           -1.265999          -1.1124704

                  resp:               resp:               resp:               resp:               resp:               resp:
        xlnhltmoney_o3            chair_o3                  o3      xlnhltmoney_o2            chair_o2                  o2
y1           .02731165            .8866291           2.0537596           1.0830145            .2473306           .77791938

                  resp:               resp:               resp:               resp:               lns1:             id1_1l:
        xlnhltmoney_o1      xlnallmoney_o1            chair_o1                  o1               _cons                  o2
y1           .05981071          -.60305853           2.4310625          -.83585784          -1.0737716           1.9340857

                id1_1l:             id1_1l:             id1_1l:              id1_1:            mnu2_1l:            mnu2_1l:
                    o1                  w1                  w2                  o3                  o2                  o1
y1          -.40709226            .6743125            .0667051            1.086956          -2.4557253           1.0983835

               mnu2_1l:            mnu2_1l:             mnu2_1:
                    h1                  h2                  o3
y1            1.657076           1.3790065           .44428856

*/









